Seasonal and geographical differences in the ruminal microbial and chloroplast composition of sika deer (Cervus nippon) in Japan

To understand the nutritional status of culled wild sika deer (Cervus nippon), we compared the ruminal microbes of deer living in habitats differing in food composition (Nagano winter, Nagano spring, and Hokkaido winter) using next-generation sequencing. Twenty-nine sika deer were sampled. Alpha and beta diversity metrics determined via 16S and 18S rRNA amplicon-seq analysis showed compositional differences. Prevotella, Entodinium, and Piromyces were the dominant genera of bacteria, fungi and protozoa, respectively. Moreover, 66 bacterial taxa, 44 eukaryotic taxa, and 46 chloroplastic taxa were shown to differ significantly among the groups by the linear discriminant analysis effect size (LEfSe) technique. Total RNA-seq analysis yielded 397 significantly differentially expressed transcripts (q < 0.05), of which 48 (q < 0.01) were correlated with the bacterial amplicon-seq results (Pearson correlation coefficient > 0.7). The ruminal microbial composition corresponded with the presence of different plants because the amplicon-seq results indicated that chloroplast from broadleaf trees and Stramenopiles-Alveolates-Rhizaria (SAR) were enriched in Nagano, whereas chloroplast from graminoids, Firmicutes and the dominant phylum of fungi were enriched in Hokkaido. These results could be related to the severe snow conditions in Hokkaido in winter and the richness of plants with leaves and acorns in Nagano in winter and spring. The findings are useful for understanding the nutritional status of wild sika deer.


Results
Richness and evenness of ruminal microbial and chloroplast composition among sampling seasons and locations revealed using amplicon-seq. For amplicon-seq, raw read data were obtained from 29 deer [12 Nagano winter (NW) deer, 12 Nagano spring (NS) deer, and 5 Hokkaido winter (HW) deer]. Two read datasets, both from NS deer, were discarded because too few reads met the quality filtering criteria. The profiles of 27 sika deer that were enrolled and 2 that were excluded based on quality checks in this study are shown in Supplementary Tables S1 and S2. The total read count (mean read count per sample) of this study was 7,042,740 (260,842 ± 87,194 SD), 6,387,893 (236,588 ± 138,486 SD), and 2,151,717 (79,693 ± 75,712 SD) for Bacterial 16S rRNA (Bac), Eukaryal 18S rRNA without chloroplasts and host mitochondria (Euk), and Chloroplast 18S rRNA (Chl) gene sequences, respectively. The total numbers of operational taxonomic units (OTUs) (mean OTUs observed per sample) identified in this study were 10,055 (372 ± 42 SD) Bac OTUs, 1252 (46 ± 14 SD) Euk OTUs and 1006 (37 ± 12 SD) Chl OTUs.
The mean read count per sample, OTUs per sample, alpha diversity and related statistics for season and location groups (NW, NS, and HW) of the ruminal microbial community diversity including Bac, Euk, and Chl are summarized in Table 1. For all Bac, Euk, and Chl amplicon-seq data, rarefaction curves were constructed to ensure a sufficient sequencing depth (89,000, 24,000, and 9000 reads, respectively) for evaluating the dominant ruminal microbial and chloroplast composition (Fig. 1). In the Bac analysis, significant differences were not observed among the groups; however, in the Euk analysis, significant differences were observed among the groups in terms of the Shannon index (p < 0.01) and Pielou's evenness (p < 0.01). In particular, the NW group showed the lowest eukaryotic richness (q < 0.05) and evenness (q < 0.05). In the Chl analysis, significant differences among the groups were observed only for observed OTUs (p < 0.05) and the Shannon index (p < 0.05) ( Table 1).
To examine the factors affecting ruminal microbial and chloroplast richness, a generalized linear model was applied to the OTUs identified in 27 sika deer based on Bac, Euk, and Chl sequences, including groups, sex, and age as factors. The group factor (NS) significantly predicted Chl diversity (p < 0.05). The analysis of deviance also showed significant differences only for groups in the Chl analysis (p < 0.05) and not for age and sex in the Bac, Euk, and Chl analyses (Supplementary Table S3). Table 2, between NS and HW, statistically significant differences were observed in Jaccard (q < 0.01), Bray-Curtis (q < 0.01), and unweighted UniFrac (q < 0.01) distances for Bac, Jaccard (q < 0.01) and unweighted UniFrac (q < 0.05) distances Table 1. Number of read counts, OTUs, and results of alpha-diversity metrics in the ruminal microbial and chloroplast composition of sika deer summarized by season and location. a NW, Nagano winter; NS, Nagano spring; HW, Hokkaido winter. b () indicates standard deviation. c q-values for each pair were calculated by pairwise Kruskal-Wallis test. q-values in NW, NS, and HW show the statistics against HW, NW, and NS, respectively. d p-values for all groups were calculated by Kruskal-Wallis test. e Mitochondria and chloroplast taxa were excluded.  Taxon-level differences in ruminal microbial and chloroplast composition detecting using amplicon-seq. The linear discriminant analysis (LDA) effect size (LEfSe) technique was used to compare taxon abundances among sampling seasons and locations at the species level (level 7). The LEfSe technique identified 66 Bac taxa, 44 Euk taxa, and 46 Chl taxa with LDA scores above 3 and significant differences (p < 0.05) among NS, NW, and HW (Supplementary Tables S10 to S12). Relevant features of these taxonomic or phylogenetic trees were visualized by cladograms (Figs. 7,8,9).

Differences in ruminal microbial and chloroplast composition among sampling seasons and locations revealed using amplicon-seq. For the beta diversity metrics shown in
Of the 66 Bac taxa, 27 showed a significant difference (p < 0.01) among groups ( Fig. 7 and Supplementary Table S10). Of the 27 Bac taxa, the taxa identified in HW were g_Lachnospiraceae XPB1014 group, g_Lachnoclostridium 10, g_Ruminococcus 1, s_Fibrobacter UWB11, s_F. aceaebacterium, and s_F. succinogenes subsp. succinogenes, whereas those in NS were g_Ruminobacter and s_Treponema AC3 and those in NW were f_Ruminococcaceae, f_Oligoflexales 0319 6G20, and s_Ruminococcaceae UCG_014.
Of the 46 Chl taxa, 13 showed a significant difference (p < 0.01) among groups ( Fig. 9 and Supplementary Table S12). Of the 13 Chl taxa, the taxa identified in HW were o_Poales and g_Panax, whereas those in NS were o_Fabales, o_Laurales, g_Cinnamomum, and C. camphora.   Table 3, and the raw numerical data for the DE analysis are provided in Supplementary Table S13.   Red, green, and blue indicate significantly different groups (Hokkaido winter, Nagano spring, and Nagano winter, respectively), with the diagram and species classification at the phylum, class, order, family, and genus levels shown from the inside to the outside. Yellow circles represent species with no significant difference. Letters in the cladograms indicate the significantly different taxa in each clade including more than 2 taxa. Eukaryotic taxa excluded chloroplasts and host mitochondria from 18S rRNA amplicon-seq analysis. A total of 19 taxa with differential expression and 48 DE transcripts were detected to have at least 1 correlation and covariance of more than 0.7 ( Fig. 10 and Supplementary Table S14). Almost all taxa with differential expression and DE transcripts were in a comparison with HW ordered first (17 and 46, respectively) when grouped by TCC-baySeq. The abundant taxa with differential expression that correlated with the DE transcripts (number) were as follows: g_Prevotella 1 (3), g_Fibrobacter (2), and g_Lachnospiraceae NK4A136 group (2). The abundant functions and pathways (gene symbols) of the DE transcripts correlated with the taxa showing differential expression were as follows: the cellular component of flagellin (FLA, FLA2, FLAB1, and FLAB2), the Embden-Meyerhof-Parnas (EMP) glycolysis pathway (ENO, G3P, GLGC, G6PI, and PGK), pyruvate fermentation to lactate (LDH), and phosphoenolpyruvate:carbohydrate phosphotransferase (PTS) (PT1 and PTHP).
Correlation between the results of 16S rRNA amplicon-seq and PICRUSt2 analysis. The metagenomic functions of Bac taxa were predicted using PICRUSt2. A total of 383 pathways (mean 324 ± 17 SD pathways per sample) related to the MetaCyc pathway were identified in PICRUSt2 through structured mappings of Enzyme Commission (EC) gene families to pathways. The raw MetaCyc pathway data were used for the association analysis because the abundances had been normalized by the PICRUSt2 pipeline.
A total of 96 correlations (14 taxa with differential expression and 52 MetaCyc pathways) were detected to have at least 1 correlation and covariance of more than 0.7 ( Fig. 11 and Supplementary Table S15). Of the 14 taxa with differential expression, the order predicted by TCC-baySeq was HW > other (number = 6 taxa), NW > other (1 taxon), other > HW (5 taxa), and other > NW (2 taxa). Of the 96 correlations (the number of each MetaCyc class of degradation/utilization/assimilation, biosynthesis, generation of precursor metabolites and energy, and macromolecule modification), 49 were included in the taxa for HW > other (33, 16, 0, 0), 10 for other > NW (0, 10, 0, 0), 34 for other > HW (11,17,5,1), and 3 for NW > other (2, 0, 1, 0). The probabilities of MetaCyc classes between the HW > other and other > HW comparisons were significantly different according to Fisher's exact test (p < 0.01). In the HW > other comparison, aromatic compound degradation was the most abundant (the total number of correlations was 26) degradation/utilization/assimilation function. In the other > HW comparison, cofactor, carrier, and vitamin biosynthesis associated with biosynthesis of cytoplasmic membrane were the most abundant (the total number of correlations was 8) biosynthesis functions.

Discussion
We examined the ruminal microbial and chloroplast composition of wild sika deer using next-generation sequencing and subsequent analysis. Red, green, and blue indicate significantly different groups (Hokkaido winter, Nagano spring, and Nagano winter, respectively), with the diagram and species classification at the phylum, class, order, family, and genus levels shown from the inside to the outside. Yellow circles represent species with no significant difference. Letters in the cladograms indicate the significantly different taxa in each clade including more than 2 taxa. The plant taxa were obtained from 18S rRNA amplicon-seq analysis.  16,17,19 and related species (wild roe deer 13 , elk 14,18 , reindeer 20 , ruminants 6 ). Of the previous reports, Henderson et al. 6 reported the dominant core rumen bacteria of 32 species of ruminants and a total of 742 individual animals. The core bacteria were similar across all ruminants; the top mean relative abundances of bacterial groups in 60 deer were reported here to be almost 70% and were observed for Prevotella, unclassified Clostridiales, Bacteroidales, Ruminococcaceae, Lachnospiraceae, Veillonellaceae, and Ruminococcus.
In the analysis of the dominant Euk taxa (Fig. 5B and Supplementary Table S7), g_Piromyces and g_Entodinium were the dominant genera. This dominance has been reported in ruminants 6 , sika deer in Japan 10 and China 17 , and red deer in New Zealand 32 . Henderson et al. 6 reported that the abundances of the dominant core rumen protozoa in the deer were almost 40% for Entodinium and 10% for Metadinium, Polyplastron, and Epidinium. This percentage for Entodinium is similar to our result for NS and NW but not HW. Of the 60 deer samples, five samples of Japanese sika deer (Cervus nippon yesoensis) obtained postmortem in January of the winter season in Hokkaido, Japan, were included in the report of Henderson et al. 6 . Similar to our results for HW, the average number of protozoa sequence reads was less than 10, while that of bacteria was almost 10 thousand.  Table 3. Group-specific predicted taxa in the ruminal microbes annotated by Swiss-Prot and analyzed by TCC-baySeq (q-value < 0.05).

Figure 10. Pearson correlation coefficients of rumen microbial composition between amplicon-seq and total
RNA-seq data. The taxonomy_7_levels.txt file in the SILVA database was used for annotation of bacterial OTUs in the amplicon-seq data. BLAST + against the Swiss-Prot database was used for the annotation of expected ORF sequences in the total RNA-seq data. Correlations with coefficients between 0.7 and 1.0 are shown. The color of the taxa with differential expression and gene symbols indicate statistically significant differences among the groups (q < 0.01) analyzed by TCC baySeq. The flagellin cellular component and EMP glycolysis pathway are shown as filled circles and triangles, respectively. www.nature.com/scientificreports/ Thus, we considered the results of our study to be informative for understanding the diversity of ruminal microbes and chloroplasts associated with different food habits developed over a long history of adaptation to natural habitats in Japanese sika deer. In addition, they could be useful for understanding the rumen composition of livestock ruminants. Genera enriched in the rumen among the groups. The groups differed in ruminal microbe and plant composition, with an order from most to least diverse of NW, NS, and HW, as shown by alpha diversity metrics for Euk and Chl and beta diversity distance for Bac and Euk (Figs. 1 and 3, and Tables 1 and 2).
Of the Bac taxa, the LEfSe technique revealed that g_Prevotella 1 and Succinivibrionaceae were enriched in NS ( Fig. 4A and Supplementary Table S7). Prevotella ruminicola, P. albensis, P. brevis, and P. bryantii have been reported to have expanded ruminal niches under starch-, pectin-, xylan-, and sugar-rich conditions 5 . It has been reported that Prevotella and Succinivibrionaceae are more abundant in animals fed diets containing concentrate than in those fed forage-rich diets 6 . These taxa may be major producers of propionate and the propionate precursor succinate and thus might be responsible for the greater levels of propionate formed under concentraterich diets 6 . It has been reported that the proportion of Prevotella and the ruminal concentration of ammonia decreased on day 28 after the start of cellulose-rich corn silage feeding in sika deer 16 . A positive association between Prevotella and the ruminal concentration of butyrate was also shown previously 16 . On the other hand, in the present study, cellulolytic bacteria such as Lachnospiraceae and Ruminococcaceae 33 were enriched in HW ( Fig. 7 and Supplementary Table S10). In particular, R. flavefaciens and F. succinogenes, enriched in HW, are well-studied cultivable bacteria 5,34,35 . Consistent with our results, it has been reported that the sika deer rumen flora contains a higher proportion of R. flavefaciens in winter than in summer at Shiretoko 11 . Under winter conditions in the same habitat, a shortage of food corresponds to a decrease in ruminal concentrations of total VFAs, propionate, butyrate, ammonia and mineral contents in sika deer 10,36 . Given the accordance between the results of previous reports and those of the present study, it is hypothesized that the amounts of nutrients in the rumen, such as cellulose, ammonia and total VFAs, are associated with the proportion of Bac taxa enriched in Nagano and Hokkaido, such as Prevotella and cellulolytic bacteria, respectively ( Fig. 7 and Supplementary Table S10). www.nature.com/scientificreports/ Among Euk taxa, the relative abundance of Entodinium in HW was markedly low in comparison to NW and NS ( Fig. 8 and Supplementary Table S11). Our results are consistent with the fact that Entodinium density in winter decreases to 1/10 of that in summer at Shiretoko 10 . In winter, due to the limited availability of deciduous broadleaf tree leaves, sika deer eat foods such as ground plants, bark, twigs, and lichens, which are less abundant and have lower-quality crude protein than tree leaves 10 . The sika deer in HW therefore could be considered adapted to utilizing low-quality diets because rumen protozoa have been reported to make minor contributions to fiber digestion 10 .
The whole transcriptomes of Entodinium caudatum have been reported using RNA-Seq analysis 37 . E. caudatum is considered to prefer starch since it has more transcripts involved in the use of starch than in the use of cellulose and hemicellulose. In addition, E. caudatum ferments sugars to VFAs. In the present study, C. glandium was enriched in NW ( Fig. 8 and Supplementary Table S11). In winter in Nagano, sika deer might prefer to eat acorns because weevil larvae often grow in acorns over the winter 38 . Acorns themselves are rich in starch and lipids 39 . The ruminal concentration of total VFAs was reported to be higher in summer and autumn than in winter 10 . Rumen protozoa affect the meat quality of ruminant mammals in terms of fatty acid composition. It has been reported that rumen protozoan membranes contain high concentrations of conjugated linoleic acids (CLAs) and vaccenic acid 40 , and the total ciliate number in the rumen has a positive correlation with cis-9, trans-11 CLAs and vaccenic acid depositions in lamb meat 41 . These findings suggested that rumen Entodinium could play an important role in storing nutrients in the muscle for severe winter seasons after deer consume fallen leaves and acorns in autumn.
Among the Euk taxa, the anaerobic fungi f_Neocallimastigaceae and g_Piromyces, well-established species with utility for lignocellulose bioprocessing 42 , were more enriched in HW than in NW and NS ( Fig. 8 and Supplementary Table S11). It has been reported that the transcripts of biomass-degrading enzymes are repressively controlled by bacteria via glucose stimulation 42 . The comparison of transcriptomes among anaerobic and aerobic fungi and rumen and nonrumen bacteria showed that P. rhizinflata had almost twice as many putative pectinases as the other anaerobic lignocellulolytic fungi 43 . These putative pectate lyases were conserved among the anaerobic fungi, whereas these enzymes were much less commonly found in the other groups of organisms, particularly bacteria 43 . Under winter conditions at Hokkaido, nutritional shortages could be required for Piromyces to degrade pectin within the bark. In our results, putative glycolysis-and glycogenesis-estimated DE transcripts were observed in the rumen bacteria of HW (Table 3 and Fig. 10). It was considered that the number and function of rumen lignocellulolytic fungi might be regulated by cellulolytic bacteria to overcome severe winter conditions. Seasonal differences in the food habits of wild sika deer in Japan. The alpha and beta diversity metrics of Chl taxa were different among NW, NS, and HW, suggesting seasonal changes in the food habits of sika deer in Japan (Figs. 1, 2, 3 and Tables 1 and 2). The present study indicates that the NW diet is enriched in broadleaf trees (Fagales, Morus, and Ilex), which bear many acorns and winter buds ( Fig. 6 and Supplementary  Table S9), whereas the NS diet is enriched in spring-to summer-flowering deciduous broadleaf trees (Populus, R. pseudoacacia, and A. quinate), evergreen trees (C. camphora, Fokienia) and forbs (P. tenella) (Fig. 6 and Supplementary Table S9). The food habits of sika deer living close to Nagano have been described 24,25 . In winter to spring in the wildlife reserve area of Ashio, Tochigi, the predominant rumen contents of sika deer were culm sheaths, dead leaves, barks and twigs, conifers, and graminoids (from most to least prevalent) 24 . In autumn to winter at Mount Akagi, Gunma, the diet was composed of graminoids, dead leaves of broadleaf trees, nuts, and berries 25 . These results indicated that the sika deer in Nagano consume mostly broadleaf tree leaves and nuts during winter to spring.
In the present study, the rumen composition of HW was enriched in graminoids (Poales) and forbs (Panax) ( Fig. 6 and Supplementary Table S9). In winter, at Shiretoko, Hokkaido, bark and twigs, graminoids, and leaves were the main foods in the rumen contents of sika deer 10 . In February, at Ashoro and Onbetsu, Hokkaido, the rumen contents were composed of broadleaf tree leaves, bark and twigs, graminoids and forbs 23 . Among the Poales, Sasa bamboo (Sasa nipponica) has been reported to be a predominant food species, especially in winter, for wild sika deer in Japan 10,23-25 . S. nipponica has evergreen leaves, and its nutritional content is moderate compared with that of other plants 23 . The Chl taxa enriched in HW indicated that Poales taxa are an important source of nutrients for sika deer in Hokkaido during severe winter periods with heavy snow.
The associations among nutritional sources, fermentation, and rumen microbial and chloroplast composition. Given the discussion above, it is hypothesized that the rumen bacteria (e.g., Lachnospiraceae and Ruminococcaceae) and anaerobic fungi (P. finnis) enriched in HW prefer to grow in cellulose-rich and starch-and lipid-poor conditions, while abundant starch and lipids are fermented by Prevotella, Succinivibrionaceae and ciliates (e.g., Entodinium) to generate large energy stores. It is suggested that the rumen bacteria and anaerobic fungi in HW extract glucose and xylose from the fermentation of cellulose, hemicellulose and lignin while saving and using small energy sources. It has been reported that the number of rumen bacteria increases and that of Entodinium decreases in winter based on morphological examinations of deer 10,44 . Significant numbers of ruminal bacteria can be consumed by protozoa, resulting in an inverse relationship between protozoal and bacterial densities 5 . Total RNA-seq analysis revealed many glycolysis-and glycogenesis-estimated DE transcripts in the rumen in the HW > other comparison (Table 3). Furthermore, the majority of the taxa with differential expression that were highly associated with the predicted DE transcripts were also enriched in HW (Fig. 10). The MetaCyc pathways associated with aromatic compound degradation were correlated with the bacterial taxa identified in the HW > other category (Fig. 11). Aromatic compounds are secondary metabolites of lignin, so this result suggested that sika deer consume low-nutrient foods in HW 45  www.nature.com/scientificreports/ bacterial activity was enhanced under rumen conditions in HW. The decrease in protozoa indicated that the interaction of rumen bacteria with anaerobic fungi compensates for the presence of protozoa.
Limitations. We could not examine fermentation parameters in the rumen or the quality of meat. Nagano sika deer and Hokkaido sika deer are different subspecies (C. n. centralis and yessoensis, respectively), and an effect of host genetic background on the rumen microbiome and metabolites in sika deer and elk hybrids has been reported 15 . Further research is warranted to examine the roles that differences in the season and location play in structuring the rumen microbiome in the same genetic background, determining rumen physiological status and driving the quality of game meat from sika deer in Japan.
To reduce contamination, we employed RNAlater to stabilize the samples and to prevent bacterial multiplication 46 . However, environmental contamination could potentially have been present due to the lack of a negative control in our study. Although a negative control would not completely ensure the absence of possible environmental contamination in the samples 47 , bioinformatic tools such as DECONTAM would be helpful for identifying and removing contaminant-related DNA sequences in microbiome data via statistical classification procedures 47 . Additional studies using negative controls and DECONTAM could improve the accuracy of the profiling of taxa with low abundance.
Correlation between total RNA-seq and amplicon-seq data. Strong correlations between the annotated DE transcripts and bacterial taxa with differential expression were observed in the present study (Fig. 10, and Supplementary Table S14). The correlated taxa with differential expression were consistent with the genera detected as enriched by the LEfSe technique, such as f_Lachnospiraceae (E. ruminantium), g_Lachnoclostridium 10, g_Lachnospiraceae XPB1014 group, g_Fibrobacter, and g_Ruminococcus 1 in HW and g_Prevotella 1 in NS (Supplementary Tables S7 and S14). The correlated flagellin genes could indicate that rumen bacteria have motility and can adhere to plant tissues (Fig. 10) 48 . The correlation with the EMP pathway, LDH, PTS, and sucrose catabolism could indicate that rumen bacteria utilize glucose, sucrose, and pyruvic acid for glycolysis, resulting in the production of VFAs ( Fig. 10 and Supplementary Table S14) 49,50 . Biological information converted from annotated transcripts has been used to study the seasonal prevalence of functional marine microorganisms 27 . The large set of gene expression data from the total RNA-seq analysis could indicate the accuracy of genomic data from amplicon-seq from different perspectives. One of the criticisms of total RNA-seq is that it is limited in the biological information provided by annotated transcripts. In this study, only 119 of the 397 DE transcripts provided information by a homology search using BLAST + against the Swiss-Prot database. This limitation is decreasing over time as the amount of biological information in the database continues to grow.
In this study, MetaCyc pathway abundances revealed by PICRUSt2 analysis were used to analyze correlations with amplicon-seq data 51 . PICRUSt2 predicts the functional potential of the bacterial community using 16S rRNA marker gene data. One of the limitations of this amplicon-based approach is that PICTRUSt2 cannot provide sufficient resolution to distinguish strain-specific functionality 51 . In this study, the biological function associated with amplicon-seq results differed between the total RNA-seq and PICRUSt2 approaches. Total RNA-seq may be superior to amplicon-seq when information is required from a live subject, such as the detection of the EMP glycolysis pathway and locomotion for rumen bacteria. Moreover, it has been reported that both approaches show varying results 30,52 , so a combination of different approaches is desirable to permit an evaluation of the accuracy of rRNA data.

Conclusion
In the present study, using 16S and 18S rRNA amplicon-seq and total RNA-seq analysis, we revealed differences in population-level microbiome and transcriptome diversity between seasons (spring/winter) and locations (Nagano/Hokkaido) in wild sika deer in Japan. Diversity was reflected by alpha diversity metrics for Euk and beta diversity metrics for Bac and Euk. The LEfSe results showed the enrichment of specific Bac and Euk taxa, and a difference in food habits was revealed by the enriched Chl taxa. Previous reports and the results of the present study provide insight into the relationship between the ruminal community (bacteria, fungi, and protozoa) and food habits (limited vs. abundant feeding conditions). The strong correlation between the results of amplicon-seq and total RNA-seq analysis showed the accuracy of the present data and their ability to provide information about the live cell function (locomotion ability and glycolysis) of ruminal bacteria. These data are useful not only for understanding the nutritional status of wild sika deer in Japan but also for understanding microbial community dynamics during the fermentation process in ruminants. Tables S1 and S2) , which were close to the sampling area. In Hokkaido, hunted sika deer were immediately dressed, and their rumen samples were collected in the field. In Nagano, the culled deer were moved to a slaughter facility operated by the Komoro city government. We used RNAlater (Thermo  -GTC TCG TGG GCT  CGG AGA TGT GTAT AAG AGA CAG GAC TAC HVG GGT ATC TAA TCC-3′) 53 . The first forward PCR primer targeting 18S rRNA (V4 region) was modified TAReuk454FWD1 (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCA GCA SCY GCG GTA ATT CC-3′), and the reverse primer was TAReukREV3 (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GAC TTT CGT TCT TGA TYR A-3′) 54 . The first PCR products were purified using Agencourt AMPure XP (Beckman Coulter, Brea, CA, USA) and eluted with 50 µl of water. A second PCR was performed using the primers of a Nextera XT index kit (Illumina). The cycle conditions were 98 °C for 2 min; 12 cycles at 98 °C for 15 s, 55 °C for 15 s, and 68 °C for 1 min; and then 68 °C for 6 min. The concentration of PCR products was measured by a Qubit DNA HS Kit (Thermo Fisher) and equalized by water.

Animals and sampling (Supplementary
A total RNA-seq library was constructed using a SMARTer stranded RNA-seq kit (Clontech, San Francisco, CA, USA). PCR amplification was repeated for 12 cycles according to the manufacturer's instructions. PCR products were eluted with 15 µl of water, measured with a Qubit DNA HS Kit and equalized with water.
Quantification of both pooled sequencing libraries was performed by qPCR using the Library Quantification Kit for Illumina sequencing platform (KAPA Biosystems, Boston, USA). The concentration was adjusted to 2.0 nM. Paired-end reads 250 base pairs (bp) in length with 5% PhiX spike were loaded for Illumina MiSeq sequencing.
Sequence analysis. The obtained 16S and 18S rRNA gene sequence reads were analyzed using the Quantitative Insights Into Microbial Ecology (QIIME2) pipeline 55 with DADA2 56 for quality control. The generated OTUs were assigned using a naïve Bayes classifier for annotation 57 . Each target region was specified by primer sequences to strains in the naïve classifier with silva132_99.fna of the Silva database, release 132 58 . Annotation was performed on taxonomy_7_levels.txt in the same database. When the identified plant species were not distributed in Japan, genus-level taxa were reported in the manuscript (raw data are shown in Supplementary  Table S9).
Total RNA-seq analysis was performed according to the previously reported ARI-seq analysis pipeline 59 . In detail, the obtained total RNA-seq reads were trimmed by Trimmomatic-0.39 60 with the option "ILLUMINA-CLIP: TruSeq_LT_HT.fa:5:30:7 MINLEN:100 HEADCROP:6 LEADING:20 TRAILING:20". PhiX sequences were removed by the USEARCH 11.0.667-filter_phix option 61 . Low-complexity filtering was performed with the USEARCH 11.0.667 -filter_lowc option. The cleaned reads were assembled by Trinity v2.11.0 with a mini-mum_contig_length of 500 62 . The reads were sorted into rRNA and non-rRNA reads using SortMeRNA 63 with the reference sequences of silva-arc-16s-id95.fasta, silva-bac-16s-id90.fasta, and silva-euk-18s-id95.fasta. The sorted rRNA reads were mapped against Trinity output (Trinity.fasta) by Bowtie2 v.2.3.5.1-linux-x86_64 64 with options-U and -local. The resulting SAM files were transformed into bam files using SAMtools and then sorted. The sorted BAM files were used to calculate counts by eXpress v.1.5.1-linux_x86_64 65 . The count data were truncated with a custom script to remove reads with fewer than 10 observations. The expressed genes were annotated for the expected ORF sequence by TransDecoder software. The resulting expected ORF sequences were annotated by BLAST + against the Swiss-Prot database.
Statistical analysis. In amplicon-seq analysis, based on the taxonomy generated for microbes, we excluded reads originating from archaea, chloroplasts, and host mitochondria. The reads from the chloroplasts were used Scientific Reports | (2022) 12:6356 | https://doi.org/10.1038/s41598-022-09855-w www.nature.com/scientificreports/ for ruminal plant taxa. We removed singleton features and elements with fewer than 10 reads. The filtered feature tables and computed relative abundance per taxonomic level values were used for diversity analyses. The alpha and beta diversity metrics and their group significance were calculated using the QIIME 2 pipeline 55 . Rarefaction-curve analyses of the obtained data were performed to estimate the completeness of taxonomic community sampling. Subsequently, alpha diversity metrics (Shannon's diversity index, observed OTUs, Faith's phylogenetic diversity, and Pielou's evenness) were examined to detect differences in the number of species and species richness per sample, and statistical significance among the groups was computed using the Kruskal-Wallis (KW) test. The effects of the other variables (sex and age) among groups were examined using generalized linear models with R functions 66 . Beta diversity metrics (Jaccard, Bray-Curtis, unweighted UniFrac, and weighted UniFrac distances) were visualized using principal coordinate analysis (PCoA) with Emperor, and significant differences among and within the groups were detected using ANOSIM. The Benjamini-Hochberg false discovery rate (FDR) correction was then used to identify each metric differently represented between groups. p-and q-values < 0.05 were considered statistically significant. The LEfSe approach 67 was implemented to identify the microbial taxa that were significantly associated with the groups. The LEfSe algorithm consisted of a KW rank sum test for detecting differences between classes and LDA for detecting differences in the relevant features. The parameters were set at p = 0.05 and LDA score = 3.0 for the computation. All components of the total RNA-seq analysis were performed using R functions 66 . Tag Count Comparison (TCC) baySeq 68 was used for normalization and differential expression analysis of the multigroup RNA-seq count data. The TCC package was generated from original TbT methods (TMM-baySeq-TMM pipeline), consisting of a combination of the trimmed mean of M values (TMM) normalization 69 in edgeR 70 and DE gene detection in baySeq 71 . In this strategy, normalization of count data and DEG estimation are iterated to avoid false positives; we repeated this cycle three times in this study 72 . Transcripts were considered DE between groups at a FDR (q-value) lower than 0.05.
PICRUSt2 analysis was performed with dada2 output during QIIME2 analysis 51 . Representative sequences and feature tables were transformed by QIIME2 and then used for PICRUSt2 analysis. The default full pipeline command "picrust2-pipeline.py" was used with default parameters for this analysis. A sequence table normalized by predicted 16S copy number abundances was used as a predicted MetaCyc pathway abundance table for further correlation analysis.
The Pearson correlation coefficients between amplicon-seq data, transcriptome data and MetaCyc pathway data were calculated using the "cor" function of R 66 . The obtained matrix was used to draw a heatmap with the "gplots" package of R. In this heatmap, we visualized correlations with a coefficient between 0.7 and 1.0.

Data availability
The datasets generated and/or analyzed during the current study are available in the DDBJ Sequence Read Archive repository, accession number DRA011262 for 16S the rRNA gene sequence library, DRA011263 for the 18S rRNA gene sequence library, and DRA011264 for the total RNA-seq library.